## Replication for Unilateral Powers, Public Opinion, and the Presidency by
## Andrew Reeves and Jon Rogowski Use this code to replicate
## Table 4 and Figure 2 in the article and Table 4 of the Online SI.
rm(list = ls())
library(Zelig)
library(foreign)
library(survey)
library(plyr)
library(questionr)
library(car)
library(psych)
library(apsrtable)



## Agree = support of president
recodeTapsPos <- function(x) {
  tmp <- rep(NA, length(x))
  tmp[x == "strongly agree"] <- 4
  tmp[x == "somewhat agree"] <- 3
  tmp[x == "somewhat disagree"] <- 2
  tmp[x == "strongly disagree"] <- 1
  return(tmp)
}
## Disagree == support of president
recodeTapsNeg <- function(x) {
  tmp <- rep(NA, length(x))
  tmp[x == "strongly agree"] <- 1
  tmp[x == "somewhat agree"] <- 2
  tmp[x == "somewhat disagree"] <- 3
  tmp[x == "strongly disagree"] <- 4
  return(tmp)
}
flipScale <- function(x) {
  tmp <- as.numeric(x)
  tmp <- recode(tmp, "1=5;2=4;3=3;4=2;5=1")
  return(tmp)
}

#ssiData <-
#  ssiData[,c(
#    "demGender", "demPidNoLn", "demPidLean", "demPidClos", "demInc", "demRace", "nr2", "demPolIdeo", "demEduFull", "rog1", "rog2", "rog3", "rog4", "rog5", "rog6",
#    "rog7", "rog8", "rog9", "wts"
#  )]

load("rr_ssi.Rdata")
## Rule of Law
ssiData$law.scale <-
  alpha(data.frame(
    ssiData$rog5, ssiData$rog6, ssiData$rog7, ssiData$rog8, ssiData$rog9
  ))$score

detach("package:psych", unload = TRUE)

## Gender
ssiData$male <- rep(0, nrow(ssiData))
ssiData$male[ssiData$demGender == 1] <- 1
#ssiData$male[ssiData$gendersp == "Refused"] <- NA

## Party
ssiData$pid <- rep(NA, nrow(ssiData))
ssiData$pid[ssiData$demPidNoLn == 1 & ssiData$demPidLean == 1] <- 1
ssiData$pid[ssiData$demPidNoLn == 1 & ssiData$demPidLean == 2] <- 2
ssiData$pid[ssiData$demPidClos == 2] <- 3
ssiData$pid[ssiData$demPidClos == 3] <- 4
ssiData$pid[ssiData$demPidClos == 1] <- 5
ssiData$pid[ssiData$demPidNoLn == 2 & ssiData$demPidLean == 2] <- 6
ssiData$pid[ssiData$demPidNoLn == 2 & ssiData$demPidLean == 1] <- 7

## Income
ssiData$income <- as.numeric(ssiData$demInc)

## White
ssiData$white <- rep(0, nrow(ssiData))
ssiData$white[ssiData$demRace == 4] <- 1
## Black
ssiData$black <- rep(0, nrow(ssiData))
ssiData$black[ssiData$demRace == 3] <- 1
## Latino
ssiData$hispanic <- rep(0, nrow(ssiData))
ssiData$hispanic[ssiData$demHisp == 1] <- 1

## Obama Approval (September)
ssiData$approval <-
  recode(as.numeric(ssiData$nr2), "1=4;2=3;3=2;4=1;5=NA")

## Ideology (cons to liberal)
ssiData$ideology <-
  recode(as.numeric(ssiData$demPolIdeo), "1=7;2=6;3=5;4=4;5=3;6=2;7=1;8=NA")

## Education
ssiData$education <- as.numeric(ssiData$demEduFull)
ssiData$education <-
  recode(ssiData$education, "1=1;2=1;3=2;4=4;5=3;6=4;7=5;8=6;9=6")

## Exec Action
ssiData$execactionCong.sep.bin <-
  recode(ssiData$rog2, "1=1;2=1;3=0;4=0;5=NA")
ssiData$execactionNatSec.sep.bin <-
  recode(ssiData$rog3, "1=1;2=1;3=0;4=0;5=NA")
ssiData$execactionFeels.sep.bin <-
  recode(ssiData$rog4, "1=1;2=1;3=0;4=0;5=NA")

## Treatment Indicators
ssiData$execaction.treatment <- rep(NA, nrow(ssiData))
ssiData$execaction.treatment[!is.na(ssiData$rog1)] <- "control"
ssiData$execaction.treatment[!is.na(ssiData$rog2)] <- "congress"
ssiData$execaction.treatment[!is.na(ssiData$rog3)] <-
  "national security"
ssiData$execaction.treatment[!is.na(ssiData$rog4)] <- "pres feel"
ssiData$execaction.treatment <-
  factor(
    ssiData$execaction.treatment, levels = c("control", "pres feel", "congress", "national security")
  )



ssiData$execaction.sep <-
  rowSums(ssiData[,c("rog1", "rog2", "rog3", "rog4")], na.rm = TRUE)
ssiData$execaction.sep[ssiData$execaction.sep == 0] <- NA

ssiData$execaction.bin <-
  recode(ssiData$execaction.sep, "1=1;2=1;3=0;4=0;5=NA")
ssiData$execaction <-
  recode(ssiData$execaction.sep, "1=4;2=3;3=2;4=1;5=NA")

##SSI Models
ssiData$execaction.treatment <-
  factor(
    ssiData$execaction.treatment, levels = c("control", "pres feel", "congress", "national security")
  )
## Table 4 from main text
SSIint <-
  glm(
    execaction.bin ~  execaction.treatment  +  scale(law.scale, center = TRUE) * execaction.treatment + approval + ideology,
    data = ssiData,
    weights = wts,
    family = binomial(link = "logit")
  )
SSItable <- apsrtable(
  SSIint,
  coef.names = c(
    "Intercept",
    "President feels strongly",
    "Congress won't act",
    "National Security at issue",
    "Belief in rule of law (centered)",
    "Presidential Approval",
    "Ideology",
    "President feels strongly $\\times$ Belief in rule of law ",
    "Congress won't act $\\times$ Belief in rule of law",
    "National Security at issue $\\times$ Belief in rule of law"
  ),
  caption = "A model of support for unilateral executive action from a survey experiment manipulation different contexts.  For ease of interpretation, we standardize the measure of attitudes toward the rule of law by centering it.",
  caption.position = "below",
  label = "t:treatments",
  Sweave = TRUE,
  model.names = c("")
)

## Table 4 from Online SI
SSIint.app <-
  glm(
    execaction.bin ~  execaction.treatment +  scale(law.scale) + scale(approval)  + execaction.treatment * scale(approval) + ideology,
    data = ssiData,
    weights = wts,
    family = binomial(link = "logit")
  )
SSItable.app <- apsrtable(
  SSIint.app,
  coef.names = c(
    "Intercept",
    "President feels strongly",
    "Congress won't act",
    "National Security at issue",
    "Belief in rule of law (centered)",
    "Presidential Approval (centered)",
    "Ideology",
    "President feels strongly $\\times$ Approval ",
    "Congress won't act $\\times$ Approval",
    "National Security at issue $\\times$ Approval"
  ),
  Sweave = TRUE,
  model.names = c("")
)

## Figure 2 from the article
## Zelig is tempermental so we include only the relevant variables and drop incomplete cases.
detach("package:psych", unload = TRUE)
ssiTmp <-
  ssiData[,c(
    "execaction.bin", "ideology", "approval", "law.scale",
    "execaction.treatment", "wts"
  )]
ssiTmp$law.scale <- scale(ssiTmp$law.scale, center = TRUE)
ssiTmp <- subset(ssiTmp, complete.cases(ssiTmp))
tmp <- as.data.frame(dummify(ssiTmp$execaction.treatment))
names(tmp) <- paste("tr", gsub(" ", "\\.", names(tmp)), sep = "")
ssiTmp <- cbind(ssiTmp, tmp)
rm(tmp)

z.out <-
  zelig(
    execaction.bin ~ trpres.feel * law.scale + trcongress * law.scale + trnational.security * law.scale + law.scale + ideology + approval,
    model = "logit",
    data = ssiTmp,
    weights = ssiTmp$wts
  )
x.control <-
  setx(
    z.out, trcongress = 0, trnational.security = 0, trpres.feel = 0
  )
x.trcongress <-
  setx(
    z.out, trcongress = 1, trnational.security = 0, trpres.feel = 0
  )
x.trnatsec <-
  setx(
    z.out, trcongress = 0, trnational.security = 1, trpres.feel = 0
  )
x.trpresfeel <-
  setx(
    z.out, trcongress = 0, trnational.security = 0, trpres.feel = 1
  )

set.seed(99898)
s.control <- sim(z.out, x = x.control)
set.seed(99898)
s.congress <- sim(z.out, x = x.trcongress)
set.seed(99898)
s.natsec <- sim(z.out, x = x.trnatsec)
set.seed(99898)
s.presfeel <- sim(z.out, x = x.trpresfeel)

plot(
  rnorm(100), type = "n", ylim = c(0,.65), xlim = c(0, 3.75),
  xaxt = "n",
  ylab = "Probability of Supporting Unilateral Action",
  xlab = "Experimental Condition",
  xaxs = "i",
  yaxs = "i"
)
points(.25, mean(s.control$qi[[1]]), pch = 18)
points(1.25, mean(s.presfeel$qi[[1]]), pch = 18)
points(2.25, mean(s.congress$qi[[1]]), pch = 18)
points(3.25, mean(s.natsec$qi[[1]]), pch = 18)
segments(.25, quantile(s.control$qi[[1]], .975), .25, quantile(s.control$qi[[1]], .025))
segments(1.25, quantile(s.presfeel$qi[[1]], .975), 1.25, quantile(s.presfeel$qi[[1]], .025))
segments(2.25, quantile(s.congress$qi[[1]], .975), 2.25, quantile(s.congress$qi[[1]], .025))
segments(3.25, quantile(s.natsec$qi[[1]], .975), 3.25, quantile(s.natsec$qi[[1]], .025))
axis(
  1, at = c(.25, 1.25, 2.25, 3.25), labels = c(
    "Control", "Pres Feels Strongly", "Congress Won't Act", "Nat'l Sec at Issue"
  )
)
